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Abstract. Triangulation of a three-dimensional point from n > 2 two- 
dimensional images can be formulated as a quadratically constrained 
quadratic program. We propose an algorithm to extract candidate solu- 
tions to this problem from its semidefinite programming relaxations. We 
then describe a sufficient condition and a polynomial time test for cer- 
tifying when such a solution is optimal. This test has no false positives. 
^ ■ Experiments indicate that false negatives are rare, and the algorithm has 

excellent performance in practice. We explain this phenomenon in terms 
of the geometry of the triangulation problem. 

1 Introduction 

We consider the problem of triangulating a point X £ R 3 from n > 2 noisy 
image projections. This is a fundamental problem in multi-view geometry and 
is a crucial subroutine in all structure- from-motion systems [1] . 

Formally, let the point X £ R 3 be visible in n > 2 images. Also let Pi £ 
be a projective camera and x, £ R 2 be the projection of X in image i, i.c, 



p3x4 



Xi = nPiX,Vi = l,...,n, (1) 

where, using MATLAB notation, X = [X; l] and 77 [u; v; w~\ = [u/w; v/w\ . 

Given the set {a;, } of noise free projections, it is easy to determine X using a 
linear algorithm based on singular value decomposition (SVD) [3]. However, in 
practice we are given Xi = Xi + r/i, where rji is noise, and there are no guarantees 
on the quality of the solution returned by the linear algorithm. 

For simplicity, we assume that rji ~ A/"(0, al). Then the triangulation problem 
is to find the maximum likelihood estimate of X given the noisy observations 
{xi}. Assuming such a point X always exists, this is equivalent to solving: 

n 

argmin^ 11777^1 -£ 4 || 2 . (2) 

i 

Here and in the rest of the paper we will ignore the constraint that the point X 
has positive depth in each image. The above optimization problem is an instance 
of fractional programming which is in general hard [1] . An efficient and optimal 
solution of © is the subject of this paper. 

For n — 2, Hartley & Sturm showed that ^ can be solved optimally in 
polynomial time [3]. For n = 3, a Grobner basis based algorithm for @ was 
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proposed in [3D]. This algorithm relies on the observation that generically, the 
optimal solution to ^ is one among the 47 solutions to a certain system of 
polynomial equations. This Grobner basis method is not usefully extendable to 
higher n and efficient optimal triangulation for n > 4 views has remained an un- 
solved problem. Other approaches either use the linear SVD based algorithm as 
initialization followed by non-linear refinement which lead to locally optimal so- 
lutions with no guarantees on the run time complexity [][] , or optimal algorithms 
whose worst case complexity is exponential in n |7l8ll3j . 

Wc present a new triangulation algorithm for n > 2 views. Based on semidcfi- 
nite programming, the algorithm in polynomial time either determines a globally 
optimal solution to ^ or informs the user that it is unable to do so. Theoreti- 
cally, the operating range (in terms of image noise) of the algorithm is limited 
and depends on the particular configuration of the cameras. In practice our 
method computes the optimal solution in the vast majority of test cases. In the 
rare case that optimality cannot be certified, the algorithm returns a solution 
which can be used as an initializer for nonlinear least squares iteration. 

The paper is organized as follows. In Section[5]we formulate triangulation as a 
constrained quadratic optimization problem. We present semidefinite relaxations 
to this problem in Scction[3] We propose our triangulation algorithm in Section[¥] 
and analyze its performance. We also provide theoretical explanation for why 
the algorithm works. Section [5] presents experiments on synthetic and real data 
and we conclude in Section [BJ with a discussion. 

Notation. We will use MATLAB notation to manipulate matrices and vectors, 
e.g., A[l : 2, 2 : 3] refers to a 2 x 2 submatrix of A. P = {Pi, . . . , P n } denotes 
the set of cameras, x = [x\;...; x n ] denotes a vector of image points, one in each 
camera, and x = [x±; . . . ; x n ] denotes the vector of image observations. Both 
x and x lie in R 2n . If y 6 R m is a vector, then y = \y; l] is the homogenized 
version of y. The inner product space of k x k real symmetric matrices is denoted 
S k with the inner product (A, B) = YjKijKk^ijBij. The set S'l C S h denotes 
the closed convex cone of positive semidefinite matrices. Wc write A >- (resp. 
A >z 0) to mean that A is positive definite (resp. positive semidefinite). 

2 Triangulation As Polynomial Optimization 

With an eye towards the future, let us re-state the triangulation problem ([2]) as 
the constrained optimization problem 

argmin ^""^ \\xj — Xj || 2 , s.t.Xi = BtPiX, V£ = 1, ...,n. (3) 

In this formulation, the constraints state that each x% is the projection of X in 
image i. Let us now denote the feasible region for this optimization problem by 

V P = \x € K 2 " 3X e M. 3 s.t. x t = nP.X, Vi = 1, . . . , n} . (4) 
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For any x £ Vp, we can recover the corresponding X using the SVD based 
algorithm [3j . Now we can re-state purely in terms of x 

argmin \\x — x\\ 2 , s.t. x £ Vp. (5) 

X 

Let Vp 2 Vp be the closure of Vp, meaning that Vp contains all the limit 
points of Vp . Then consider the following optimization problem: 

argmin ||x — i|| 2 , s.t. x £ Vp. (6) 

X 

The objective function in these two optimization problems is the squared 
distance from x to the sets Vp and Vp respectively. Since Vp is the topological 
closure of Vp it can be shown that any solution x* to which is not in Vp is 
arbitrarily close to a point in Vp, and the optimal objective function values for ([5]) 
and (J6j) are the same. Thus, solving ([6]) is essentially equivalent to solving ([5]). 

The set Vp is a quasi-projective variety. A variety is the zero set of a finite 
set of polynomials, and a quasi-projective variety is the set difference of two 
varieties. Therefore, Vp is also a variety [TS]. Hey den & Astrom [5] show that 

Vp = jx £ R 2 ' 

Here fij(x) = xj FijXj = are the bilinear/quadratic cpipolar constraints, where 
Fij 6 R 3x3 is the fundamental matrix for images i and j. The second set of 
constraints tijk(xi,Xj,Xk) = are the the trilinear/cubic constraints defined by 
the trifocal tensor on images i,j and fc. 

At the risk of a mild abuse of notation, we will also use to denote a 
(2n + 1) x (2n + 1) matrix such that fij(x) = i T F^x. The construction of this 
matrix involves embedding two copies of the 3x3 fundamental matrix (with 
suitable reordering of the entries) in an all zero (2n + 1) x (2n + 1) matrix. 

Now, let Wp be the quadratic variety 

W P = {x e R 2 ™ | x T FijX = 0, 1 < i < j < n} . (8) 

For n = 2, since there are no trilinear constraints Wp = V p but for n > 3, in 
general Wp 2 Vp. For n > 4, Hey den & Astrom show that if the camera centers 
are not co-planar then Wp = Vp [3]. Note that for n = 3, the camera centers are 
always co-planar. Therefore, for n = 2, and for n > 4 when the camera centers 
are non-co-planar, we can just optimize over the quadratic variety Wp: 

argmin ||a; — x\\ 2 s.t. x £ Wp. (9) 

X 

However, we cannot just ignore the co-planar case as a degeneracy since it 
is a common enough occurrence, e.g., an object rotating on a turntable in front 
of a fixed camera. If all the camera centers lie on a plane irp, then solving ([5]) 
instead of ^ can result in spurious solutions, i.e. a projection vector x* for 
which there is no single point X* £ R 3 that projects to x* for each image i. This 



fij(x) =0, 1 < i <j < n 
tij.k(z) — 0, l<i<j<k<n 



(7) 
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can happen if each x* lies on the image of the plane up in image i. It is easy to 
reject such a spurious x* by checking if it satisfies the trilinear constraints. 

From here on, we will focus our attention on solving ([9]) and in Section [5] we 
will show that solving ((9]) instead of (|6|) is not a significant source of failures. 

Let us now define the polynomial 



x T 



I -x 



x T Gx, 



(10) 



where / is the 2n x 2n identity matrix. Observe that G £ Sj n+ , and the Hessian 



of g is V 2 g = 21. Similarly, let Fi- 



+ 

b R ij , where ^ eS 2n , 



p2n 



, and 



flij £ K. Then V 2 = 2Hij. We re- write (|9]) as the quadratically constrained 
quadratic program (QCQP) 



arg mm x 



s.t. x T Fiji = 0, 1 < i < j < n. 



(11) 



3 Semidefinite Relaxation 

We re- write ([TTjl as the following rank constrained semidefinite program (SDP). 

arg min (G, Y) 

s.t. (Fi^Y) =0, 1 < i < j < n, 

(E,Y) = 1, (12) 
Y e Sl n+1 , 

rank(F) = 1. 

Here E s S 2n+1 is an all zero matrix except for its bottom right entry which 
equals one. The problems (fTTj) and (fl~2"|) are equivalent: x is feasible (optimal) 
for (fTT|) if and only if Y = xx T is feasible (optimal) for ([T^jl . 

Solving rank constrained semidefinite programs is NP-hard [22]. Dropping 
the rank constraint gives the primal semidefinite program 



arg min (G, Y) 

s.t. (Fij,Y) =0, 1 < i < j < n, 
{E,Y) = 1, 
Y e 5 2 " +1 . 

The dual of this primal semidefinite program is 



(13) 



arg max p 

subject to G + J2 A,,/',, -pE^O, ( 14 ) 
Xij,p el, 1 < i < j < n. 



The primal SDP f| 13|) is also known as the first moment relaxation and its 
dual SDP ([T4")) is known as the first sum of squares relaxation. They are instances 
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of a general hierarchy of semidefinite relaxations for polynomial optimization 
problems [TU]. Problem (fTTj) is also the Lagrangian dual of (fTTj) . 

The remainder of this paper is dedicated to the possibility that solving the 
triangulation problem is equivalent to solving these semidefinite relaxations. Let 
us denote by g*, gi uad ) g m ° m ; g,s°s ; ^ e optimal solutions to the optimization prob- 
lems ([5]), (fTTj) . (H3J), (HH respectively. Then the following lemmas hold. 

Lemma 1. For n, g* > gi uad . for n = 2, or n > 4 wit/i non- co-planar 
cameras, g* = gi uad . 

Proof. The claim follows from the discussion in Section [5] after the definition of 
W P . □ 

Lemma 2. g^ uad > g mom . 

Proof. This is true because (fT3"j) is a relaxation of (fTTj) . □ 

Lemma 3. For all n, g mom = g s ° 8 . Moreover, there exist optimal Y*,\*j and 
p* that achieve these values. 

Proof. The inequality g mom > g sos follows from weak duality. Equality, and the 
existence of Y* , X*j and p* which attain the optimal values follow if we can show 
that the feasible regions of both the primal and dual problems have nonempty 
interiors [521 Theorem 3.1] (also known as Slater's constraint qualification). 

For the primal problem, let x G R 2 ™ be any feasible point for the triangulation 
problem ([9]) (such a feasible point always exists) and let D = diag(l, ...,1,0) G 
S 2n+1 . It is easy to show that Y = xx T + D is positive definite and primal 
feasible. For the dual problem, take \j = and p = — 1 and verify G + E >~ 0. 

□ 

4 The Algorithm and its Analysis 

We propose Algorithm [TJ as a method for triangulation. 

4.1 Correctness 

Theorem 1. Algorithm]]] terminates in time polynomial in n. 

Proof. The proof is based on three facts. One, the primal (|13[) and dual (|14l) 
have descriptions of size polynomial in n, the number of images. Two, SDPs 
can be solved to arbitrary precision in time which is polynomial in the size of 
their descriptions [22] , Three, the eigenvalue decomposition of a matrix can be 
computed in polynomial time. □ 

Before moving forward, we need the following definition. 

Definition 1. A triangulation problem is SDP-EXACT if ^ uad = 5 sos = g mom , 
i.e., the relaxations are tight. 
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Algorithm 1 Triangulation 

Require: Image observation vector x G R 2 " and the set of cameras P = {Pi, . . . , P n }- 

1: Solve the primal (|13|l and dual (|14[) SDPs to optimal Y* , \*j and p* . 

2: x = Y*[l : 2n, 2n + 1] (i.e., x is the last column of Y* without its last entry) 

3: Use the SVD algorithm to determine a world point X £ R 3 from x. 

4: if / + E KjHij y then 

5: if n = 2, or n > 4 and the cameras P are non-co-planar, then 
6: Return (OPTIMAL, X). 

7: end if 

8: if Xi = nPiX Vi = l, ...,nthen 

9: Return (OPTIMAL, X) 
10: end if 
11: end if 

12: Return (SUBOPTIMAL, X). 



We will first describe the conditions under which a triangulation problem is 
SDP-EXACT. We will then show that if Algorithm □ returns OPTIMAL then tri- 
angulation is SDP-EXACT, and further, the solution X returned by the algorithm 
is indeed optimal for triangulation. 

Theorem 2. Let x* be an optimal solution to the quadratic program The 
triangulation problem is SDP-EXACT if and only if there exist Xij £ R such that 

(i) Vg{x*) + AijV/ <3 -(&*) = and (ii) 7 + ^ A, , //, , h 0. (15) 

Before proving this theorem, we observe that it is not immediately useful 
from a computational perspective. Indeed, a priori verifying condition (i) re- 
quires knowledge of the optimal solution x* . However, the theorem will help us 
understand why the triangulation problem is so often SDP-EXACT in Section l4~51 

Let L(x,\ij,p) = g{x) + £) y X i3 f i3 (x) - p = i T (G + J2 - pE) x. 

Observe that V x L(x,Xij,p) = Vg(x) + ^ XijV fij(x) and V x L(x, Aj, ;, p) = 
2(1 + ^jXijHij). We require the following two simple lemmas, the proofs of 
which can be found in the Appendix. 

Lemma 4. If x* is the optimal solution to (jllj) and Xij satisfy condition (i), 
then L(x, Xij,g(x*)) — (x — x*) T (1 + ^2 ^ij^ij) ( x ~ x *)- Further, if condition 
(ii) is satisfied as well, then L(x, Xij, g(x*)) > 0, Vx £ R 2 ™. 

Lemma 5. If x T Ax + 2b T x + c > 0, Vx, then 

Proof (Theorem^). For the if direction, let Ay satisfy conditions (i) and (ii). 
Then from Lemma |4] we have L(x, Xij, g(x*)) > Vx G R 2 ™, which combined 
with Lemma [5] gives G + J^^ijFij — g{ x *)E >r 0. Therefore Xij and p = g(x*) 
are dual feasible, which in turn means that g sos > p — g(x*) = gi uad . Lemmas [5] 
and El give the reverse inequality, thus g sos = g mom = gq uad — g(x*). 



A b 

b T c 



y o. 
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For the only if direction, let p* and X*j be the optimal solution to the 
dual (HU). If the problem is S DP- EXACT we have the equality p* = g(x*) = g quad 
and from the dual feasibility of A*- and p* we have that G+Y^ Kj^H ~ P*E ^0. 
Taken together, these two facts imply that L(x, X*j,g(x*)) > 0, Vx G R 2n . 

Since L{x* , X*j,g(x*)) = 0, L(x, X*j,g(x*)) is a non-negative quadratic poly- 
nomial that vanishes at x* . Non- negativity implies condition (ii) (that the Hes- 
sian of the polynomial is positive scmidefinite) and the fact that zero is the 
minimum possible value of a non-negative polynomial implies that its gradient 
vanishes at x* which is exactly condition (i). □ 

Condition (ii) of Theorem [5] is automatically satisfied by any feasible Ajj 
for the dual (fT4)l . Hence, verifying SDP exactness using the dual optimal X*j 
reduces to checking condition (i) of Theorem [2] Checking this condition is not 
computationally practical, since it requires knowledge of the optimum x*. By 
slightly tightening condition (ii) we can bypass condition (i). 

Theorem 3. // {Y*, X* } , p*} are primal-dual optimal and I + J2 Kj^ij ^ 0' 
then rank(y*) = 1, and triangulation is SDP-EXACT. 

Proof. Notice that 1 + ^2 X*j Hij is the top left (2n) x (2n) block in the larger 
(2n+l) x (2n+l) positive semidefinite matrix G+J2 Kj-^ij —p*E. By hypothesis, 
I + X*jHij is nonsingular and thus has full rank equal to 2n, which implies 

rank (g + ^ X*j F i3 - p*E^ > 2n. (16) 

The dual and the primal SDP solutions satisfy complementary slackness, 
which means that (G + J2 Kj F ij ~ P* E , Y *) = °- In particular it implies that 

rank(G + ^ A*-^ - p*E) + rank(F*) < 2n + 1, (17) 

where we use the standard fact that whenever (A,B) = for A, B G <S+, then 
rank(A) + rank(B) < N. From and <[TTJ) we have rank(F*) < 1. Since 
(E, Y) = 1, we have Y* ^ and hence rank(Y~*) = 1. □ 

Line 4 of Algorithm [TJ uses Theorem [3] to establish that we have solved (1X1]) . 
Lines 5-10 of the algorithm are then devoted to making sure that the solution 
actually lies in Vp. Thus, Algorithm [1] is correct. 

4.2 Implications 

Theorem 4. // the image observations are noise free, i.e. there exists X* G M 3 
such that Xi = IlPiX* , V« = 1, . . . , n, then Algorithm^ returns OPTIMAL 

Proof. Setting Ay = satisfies the hypothesis for Theorems [5] and [3] □ 

Theorem 5. Two view triangulation is SDP-EXACT. 
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Proof. For n = 2 the triangulation problem involves minimizing a quadratic 
objective over a single quadratic equality constraint fi2(x) = 0. The conditions 
of Theorem [2] in this case reduce to finding A £ R satisfying 

Vg(x*) + XVf 12 (x*) = and J + Ai? 12 >r 0. (18) 

The existence of such a A follows directly from [14j Theorem 3.2]. □ 

Theorems [3] and [5] nearly imply that two- view triangulation can be solved 
in polynomial time. We say nearly because, despite Theorem [SJ it is possible 
that the matrix / + X*Hi2 is singular for the dual optimal A* (sec Appendix for 
an example). This is not a contradiction, since Theorem [3] is only a sufficient 
condition for optimality. Despite such pathologies, we shall sec in Section [5] that 
in practice Algorithm [1] usually returns OPTIMAL for two- view triangulation. 

4.3 Geometry of the Algorithm 

Recall that the optimization problem (jllj) can be interpreted as determining the 
closest point x* to x in the variety Wp. This viewpoint gives geometric intuition 
for why Algorithm Q] can be expected to perform well in practice. 

Lemma 6. Given x, and assuming appropriate regularity conditions at the op- 
timal solution x* to there exist A^ S R such that 

•'• • E^V/yOO (19) 

Proof. It follows from Lagrange multiplier theory |15| that there exist Lagrange 
multipliers A ?J G R such that Vg(a;*) + £ Ay V/y (x*) = 0. Observe that for (fTTI) 
Vg(x*) = 2{x* — x), which finishes the proof. □ 

If ||x* — x\\ is small, i.e. x is close to the variety Wp, then there must exist 
some Ay satisfying (fTT)|) such that ||A|| is small and hence I + J^^ij^ij ^ s a 
small perturbation of the positive definite identity matrix /. Since / lies in the 
interior of the positive semidefinite cone, these small perturbations also lie in 
the interior, that is I + ^ ^ijHij ^ 0. 

This, coupled with the fact that Ay are Lagrange multipliers at x* , yields 
the sufficient conditions in Theorem [5] for a triangulation problem to be S DP- 
EXACT. Thus, if the amount of noise in the observations is small, Algorithm [1] 
can be expected to recover the optimal solution to the triangulation problem. 
Since ify depends only on the cameras P, and not on x, the amount of noise 
which Algorithm [1] will tolerate depends only on Pj. We summarize this formally 
in the following theorem. 

Theorem 6. Let N(x*) = {x* + ^^fij(x*) I + £ Ay-ffy >- 0}. For any 
x G R 2ra , if Algorithm[l\returns OPTIMAL andx* is the optimal image projection 
vector then x £ N(x*). Conversely, ifxG N(x*) and x* is the closest point in 
Wp to x, then Algorithm[J\will return OPTIMAL. 

Proof. The proof is a straightforward application of Lemma [5] and Theorem [3] 

□ 
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5 Experiments 

Algorithm Q] was implemented using YALMIP [12] , SeDuMi [21] and MATLAB. These 
tools allow for easy implementation and the timings below are for completeness 
and should not be used to judge the runtime performance of the algorithm. 

Fundamental matrices and epipolar constraints are specified only up to scale 
and badly scaled fundamental matrices lead to poorly conditioned SDPs. This 
was easily fixed by dividing each Fij by its largest singular value. 

Algorithm [1] verifies optimality by checking the positive definiteness of I + 
^2 ^ijHij >~ 0, which requires that its smallest eigenvalue be greater than some 
S > 0. In all our experiments we set 6 = 0.05. 

Either Algorithm Q] returns OPTIMAL and the solution is guaranteed to be 
optimal, or it returns SUBOPTIMAL in which case we cannot say anything about 
the quality of the solution, even though it could still be optimal or at least serve 
as a good guess for non- linear iteration. Thus, we run Algorithm [T] on a number 
of synthetic and real- world datasets and report the fraction of cases in which the 
algorithm certifies optimality. 

5.1 Synthetic Data 

To test the performance of the algorithm as a function of camera configura- 
tion and image noise we first tested Algorithm [1] on three synthetic data sets. 
Following |16j . we created instances of the triangulation problem by randomly 
generating points inside the unit cube in R 3 . 

For the first experiment, a varying number of cameras (2, 3, 5 and 7) were 
placed uniformly at random on a sphere of radius 2. In the second experiment, 
the same number of cameras were uniformly distributed on a circle of radius 
2 in the xy-planc. In the third experiment they were restricted to the x-axis 
and were placed at a distance of 3, 5, 7, and 9 units. These setups result in 
image measurements with an approximate square length of 2 units. Gaussian 
noise of varying standard deviation was added to the image measurements. The 
maximum standard deviation was 0.2, which corresponds to about 10% of the 
image size. For each noise level we ran the experiment 375 times. Figures[lja), (b) 
and (c) show the fraction of test cases for which Algorithm [1] returned OPTIMAL 
as a function of the the standard deviation of perturbation noise. 

We first note that for n = 2 cameras, in all three cases, independent of camera 
geometry and noise, we are able to solve the triangulation problem optimally 
100% of the time. This experimentally validates Thcorcm[5]and provides strong 
evidence that for n = 2 there is no gap between Theorems [5] and [3] in practice. 

Observe that for cameras on a sphere (Figure [Ha)), the algorithm performs 
very well, and while the performance drops as the noise is increased, the drop is 
not significant for practical noise levels. Another interesting feature of this graph 
is that finding and certifying optimality becomes easier with increasing number 
of cameras. The reason for this increase in robustness is not clear and we plan 
to further pursue this intriguing observation in future research. 



10 Chris Aholt, Sameer Agarwal and Rekha Thomas 




(a) Cameras on a sphere. (b) Cameras on a circle. (c) Cameras on a line. 



Fig. 1. Fraction of synthetic experiments in which Algorithm [T] returned OPTIMAL 
versus the standard deviation of the noise level added to the images, (a) Cameras 
placed randomly on the sphere of radius 2. (b) Three cameras all placed on the xy- 
plane. (c) Three cameras placed along the x-axis. 

Cameras on a circle (Figure |Tfb)) is one of the degenerate cases of our algo- 
rithm where V p ^ Wp. We expect a higher rate of failure here and indeed this 
is the case. It is however worth noting that the algorithm does not immediately 
breakdown and shows a steady decline in performance as a function of noise like 
the previous experiment. Unlike cameras on a sphere, increasing the number of 
cameras does not increase the performance of the algorithm, which points to the 
non-trivial gap between V p and Wp for co-planar cameras. 

Finally let us consider the hard problem of triangulating a point when the 
camera is moving on a line. This case is hard geometrically and algorithmically 
as the cameras are co-planar, and this difficulty is reflected in the performance 
of the algorithm. In contrast to the previous experiments, we observe rapid 
degradation with increasing noise. 

5.2 Real Data 

We tested Algorithm [1] on four real- world data sets: the Model House, Corridor, 
and Dinosaur from Oxford University and the Notre Dame data set [T5]. Our 
results are summarized in Table [TJ 

Table 1. Performance of Algorithm [T] on real data. The column OPTIMAL reports the 
fraction of triangulation problems which were certified optimal. 



Data set 


# images 


# points 


OPTIMAL 


Time (sec) 


Model House 


10 


672 


1.000 


143 


Corridor 


11 


737 


0.999 


193 


Dinosaur 


36 


4983 


1.000 


960 


Notre Dame 


48 


16,288 


0.984 


7200 
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Model house is a simple data set where the camera moves laterally in front 
of the model house. Global optimality was achieved in all cases. 

Corridor is a geometrically hard sequence where most of the camera motion 
is forward, straight down a corridor. This is similar to the synthetic experiment 
where the cameras were all on the a;-axis. The algorithm returned OPTIMAL in 
all but one case. 

Dinosaur consists of images of a (fake) dinosaur rotating on a turntable in 
front of a fixed camera. Even though the camera configuration is hard for our 
algorithm (cameras are co- planar), global optimality was achieved in all cases. 

The three Oxford datasets are custom captures from the same camera. Notre 
Dame consists of images of the Notre Dame cathedral downloaded from Flickr. 
The data set comes with estimates of the radial distortion for each camera, 
which we ignore as we are only considering projective cameras in this paper. 
Algorithm [T] returned OPTIMAL in 98.36% of cases. 

It is worth noting here that the synthetic datasets were designed to test 
the limits of the algorithm as a function of noise. A standard deviation equal 
to 0.2 translates to image noise of approximately 10% of the image size. In 
practice such high levels of noise rarely occur, and when they do, they typically 
correspond to outliers. The results in Tabic [T] indicate that the algorithm has 
excellent performance in practice. 

6 Discussion 

Wc have presented a semidefinite programming algorithm for triangulation. In 
practice it usually returns a global optimum and a certificate of optimality in 
polynomial time. Of course there are downsides which must be taken into consid- 
eration. Solving SDPs is not nearly as fast as gradient descent methods. More- 
over, what happens in the rare cases that global optimality is not certifiable? 
Regardless, the lure of a polynomial time algorithm and a better understand- 
ing of the geometry of the triangulation problem is far too great to allow these 
hiccups to close the door on SDPs. 

Hartley and Sturm [3] show that two-view triangulation can be solved by 
finding the roots of a degree 6 univariate polynomial. For n — 3, [3U] gives a 
Grobner basis based algorithm for triangulation by finding all solutions to a 
polynomial system with 47 roots. These methods do not extend in a computa- 
tionally useful manner to n > 4. Our algorithm solves the triangulation problem 
for all values of n > 2 under the conditions of Theorem [3] when the camera 
centers are not co-planar. 

Similar to our setup, Kanatani et al. also frame triangulation as finding the 
closest point on a variety from noisy observations [§]• Unlike our work, they use 
both cpipolar and trifocal constraints, which makes the constraint set much more 
complicated. Further, they do not prove finite convergence or the optimality of 
their procedure. Our work answers their question of analyzing the shape of the 
variety to obtain a noise threshold for optimality guarantees. 
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Scmidefinite relaxations for optimization problems in multi-view geometry 
were first studied by Kahl & Henrion in [5]. They formulate triangulation as an 
optimization problem over quartic inequalities and study its moment relaxations. 
Kahl & Henrion observe that good solutions can be obtained from the first few 
relaxations. In our method the very first relaxation already yields high quality 
solutions. Further, the quadratic equality based formulation has nice theoretical 
properties that explain the empirical performance of our algorithm. 

Hartley & Seo proposed a method for verifying global optimality of a given 
solution to the triangulation problem [2J. The relation between their test and 
the definiteness test of Theorem [3] is a fascinating direction for future research. 

The study of QCQPs has a long history. There is a wide body of work devoted 
to verifying optimality of a given solution for various classes of QCQPs and using 
scmidefinite programming to solve them. Thus it is natural that statements 
similar to Theorems [5J and [3] for various subclasses of QCQPs have appeared 
several times in the past e.g. [6111117] . 

The astute reader will notice that even though Theorems [2J and [3J are pre- 
sented for triangulation, the proofs apply to any QCQP with equality constraints, 
and can be interpreted in terms of the Karush-Kuhn- Tucker (KKT) conditions 
for More recently, and independently from us, Zheng et al. have proved 

versions of Theorems [JJ and [3J for the more general case of inequality-constrained 
QCQPs [23] ■ Our formulation of triangulation as minimizing distance to a va- 
riety allows for a geometric interpretation of these optimality conditions, which 
in turn explains the performance of our algorithm. 

Finally, we use the general purpose SDP solver SeDuMi for our experiments. 
Much improvement can be expected from a custom SDP solver which makes use 
of the explicit symmetry and sparsity of the triangulation problem (e.g. for all n 
the linear matrix constraints in the primal scmidefinite program (|13[) are sparse 
and have rank at most five). 
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Appendix 

Failure of Theorem [3] for two- view triangulation 

Consider the cameras 
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where < a < b. This models the situation in which two cameras are placed 
on the z-axis at (a, 0,0) and (6,0,0), respectively, both facing the origin (i.e. 
viewing down the ir-axis). 

We can show that, given e > and x = (xi, X2) = ([0; e] , [e; 0]), the mini- 
mum value of the QCQP ([TT]) is e 2 and all points of the form 

x * = [Va^e^TO; m] ( 21 ) 

x* 2 =[e- M ; v/ M (£ - aOJ (22) 

where < fi < e arc optimal. 

The non-uniqueness of the minimizer implies that the verification matrix 
I+XH12 cannot be positive definite. So in this case, the hypothesis of Theorem[3] 
is not satisfied, although Theorem [5] still holds. 

The proof of these facts is a straightforward application of the KKT con- 
ditions. The minimum distance e 2 can be attained by setting one of the image 
points equal to the epipole, but not both. However, the minimum distance is 
also attained at points which correspond to neither image being the epipole. 

Proof of Lemma [4] 

Recall conditions (i) Vg(x*) + J2 AijV 'fy (x*) = 0, and (ii) I + Y, X v H v t 0. 
Let H = ^ijHij, b = J2 ^ijbij, and fj = ^ ^ijPij- Then 

L(x,\ij,g(x*)) = g{x)+'Yj\ijfij(x) - g(x*) (23) 

= x T (I + H) x + 2 (b - x) T x + x T x + 13- g(x*). (24) 

We wish to show that JM|) is equal to (x - x*) T (I + H) (x — x*). The fact that 
L(x,Xij,g(x*)) > is then immediate because I + H >z 0. Looking at the 
quadratic, linear, and constant terms separately, it suffices to prove 

(a) (I + H)x* =x-b and (b) x* T (I + H) x* = x T x + f3 - g(x*). (25) 

Indeed, (a) is just a restatement of condition (i), since \7g(x) = 2(x — x) and 
Vfij(x) = 2(HijX + bij). To prove (b), first recall that since x* is feasible for 
problem (JTTJ), f l] {x*) = for all In particular, /? = ~2b T x* ~x* T Hx* . Using 
this and part (a), one can verify (b) through straightforward manipulations. □ 

Proof of Lemma \E\ 

Suppose the matrix M = [A, b;b T ,c\ is not positive semidefinite, i.e. there 
exists y such that y T My < 0. Write y = [y';j] for some 7 £ K. If 7 = 0, 
then > y 1 My = y' T Ay' . But then we arrive at a contradiction by considering 
x = ny' for a scalar /igR, since for /i large enough, x T Ax + 2b T x + c < 0. 
Now if 7 7^ 0, setting x = y'/^f gives 

x T Ax + 2b T x + c = i (y T My) < 0, (26) 

which again contradicts the hypothesis. □ 



